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Complex analytical structure of Stokes wave for two-dimensional potential flow of the 
ideal incompressible fluid with free surface and infinite depth is analyzed. Stokes wave is 
the fully nonlinear periodic gravity wave propagating with the constant velocity. Simu¬ 
lations with the quadruple (32 digits) and variable precisions (more than 200 digits) are 
performed to find Stokes wave with high accuracy and study the Stokes wave approaching 
its limiting form with 27 t/ 3 radians angle on the crest. A conformal map is used which 
maps a free fluid surface of Stokes wave into the real line with fluid domain mapped 
into the lower complex half-plane. The Stokes wave is fully characterized by the com¬ 
plex singularities in the upper complex half-plane. These singularities are addressed by 
rational (Fade) interpolation of Stokes wave in the complex plane. Convergence of Fade 
approximation to the density of complex poles with the increase of the numerical preci¬ 
sion and subsequent increase of the number of approximating poles reveals that the only 
singularities of Stokes wave are branch points connected by branch cuts. The converging 
densities are the jumps across the branch cuts. There is one branch cut per horizontal 
spatial period A of Stokes wave. Each branch cut extends strictly vertically above the 
corresponding crest of Stokes wave up to complex infinity. The lower end of branch cut is 
the square-root branch point located at the distance Vc from the real line corresponding 
to the fluid surface in conformal variables. The increase of the scaled wave height H/\ 
from the linear limit 77/ A = 0 to the critical value Hmax /A marks the transition from the 
limit of almost linear wave to a strongly nonlinear limiting Stokes wave (also called by 
the Stokes wave of the greatest height). Here 77 is the wave height from the crest to the 
trough in physical variables. The limiting Stokes wave emerges as the singularity reaches 
the fluid surface. Tables of Fade approximation for Stokes waves of different heights are 
provided. These tables allow to recover the Stokes wave with the relative accuracy of at 
least 10“^®. The tables use from several poles for near-linear Stokes wave up to about 
hundred poles to highly nonlinear Stokes wave with Uc/A ^ 10“®. 
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1. Introduction 

Theory of spatially periodic progressive (propagating with constant velocity without 
change of the shape and amplitude) waves in two-dimensional (2D) potential flow of an 
ideal inc ompress i ble fluid with free surface in gravitation al field was founded in pioneering 


lyoii). ana many otners isee e.g. a dook; dv 
well as Baker & Xie (l201lll: Cowlev et al\ (Il999l 

snsKii 11 

: Grant 

ymi lor 

(T^: 

review or oiaer wc 
LonEuet-HiEEinsI ( 

)rK;s as 

20081: 

LonEuet-HiEEins & Fox (1977. 1978ll: Schwartd 

(I 1974 I); 

Tanveer 

(19911: WilliamsI (1981. 


198511 and references there in for more recent progress). There are two major approaches 


to analyze the Stokes wave, both originally developed by Stokes. The first approach is 
the perturbation expansion in amplitude of Stokes wave called by the Stokes expansion. 
That approach is very effective for small amplitudes but conver gjes very slowly ( or doe s 
not converge at all, depending on the formulation according to IPrennan et al. 1 1992l ll 
as the wave approaches to the maximum height H — Hmax (also called by the Stokes 
wave of the greatest height or the limiting Stokes wave). Here the height H is defined at 
the vertical distance from the crest to the trough of Stokes wave over a spatial period 
A. The second approach is to consider the limiting Stokes wave, which is the progressive 
wave with the highest nonlinearity. Using conformal mappings Stokes found that th e 
limiting Stokes wave has the sharp angle of 27t/3 radians on the crest ( Stokes 1880 6ll . 
i.e. the surface is non-smooth (has a jump of slope) at that spatial point. That corner 
singularity explains a slow convergence of Stokes exp ansion a s H — > Hmax- The global 
existence of the limiting Stokes wave was proven by iTolandl ( 1978l l however lacking a 
proof of a Stokes conjecture that the the jump of the slope at the crest is exactly 27t/3 
radians. The Stoke s conjecture was later independently proven bv IPlotnikovI ( 1982ll and 
Amick et al. (1982). 


It was StokesI (18806) who first proposed to use conformal mapping in order to address 
finite amplitude progressive waves. In this paper we consider a particular case of potential 
flow of the ideal fluid of infinite depth although more general case of fluid of arbitrary 
depth can be studied in a similar way. Assume that free surface is located at j/ = r]{x, t), 
where x is the horizontal coordinate, y is the vertical coordinate, t is the time and r]{x, t) 
is the surface elevation with respect to the zero mean level of fluid, i.e. r]{x, t)dx = 0. 
We consider the conformal map between the domain —oo < y ^ rj{x,t), —oo < x < oo 
of the complex plane z = x + iy filled by the infinite depth fluid and a lower complex 
half-plane (from now on denoted by C“) of a variable w = u + iv (see Fig. [Ij. The real 
line ?; = 0 is mapped into the free surface by z{w) being the analytic function in the lower 
half-plane of w as well as the complex fluid velocity potential n(r(;) is also analytic in C“. 
Both z(w) and n(ri;) have singularities in upper half-plane (here and further denoted by 
C+). 

The knowledge of singularities in C'*' would result in the efficient description of the so¬ 
lution in the physical variables. Examples of such type of solutions in hydrodynamic-type 
systems ar e numerous including e.g. the dynamics of free surface of id eal fluid with infi¬ 
nite de pth ( Kuznetsov et a6l[l993l 1994 Tanveej 1993 1 and finite depth (IPvachenko et al 


19966 1. dynamics of interface between two ideal fluids ( Kuznetsov et al .1 1 1993 ), ideal fluid 


pushe d through viscous fluid in a n arrow gap between two parallel plates (Hele-Shaw 
flow) (iMineev- Weinstein et al\ 20(10ll. the dynamics of the inte rface between ideal f luid 
and light viscous fluid ( Lushnikov 2004l and bubble pinch-off ( Turitsvn et a/.l[2009ll . In 
these systems the dynamics is determined by poles/branch cuts in the complex plane. 
Related systems correspond to the sponta neous appearance of curvature singularities on 
vortex sheets as obtained by Moord ( 1979l l. Me & Baker (1998) established that Moore’s 
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z=z(w,t) 

Figure 1. Schematic of a conformal map between the domain below the solid curve (left panel) 
in z = X + ly complex plane and the lower complex half-plane in in = u -|- in (right panel). Fluid 
occupies the domain below the solid curve in physical plane z = x -\-iy. The solid curve of left 
panel (corresponds to a free surface of the fluid) is mapped into the real line (another solid 
line) in right panel. One spatial period of Stokes wave is shown by solid lines in both panels 
in the reference frame moving with the velocity c. The dark circles mark the positions of the 
singularity closest to the fluid surface in both panels. 


singularities are present in axisymmetric vortex sheets. Inoeamov fc Oparin ( 20031) con¬ 
sidered cone-shaped nose of 27 t/ 3 degrees in axisymmetric flow. Ishihara & Kaneda (1994) 
and Hou & Hu (2003) extended Moore’s singularities to three-dimensional (3D) vortex 
sheets. Oppo site limit i s the glob al existen ce of water waves for small enough data shown 
both for 2D ( Wull200^ and 3D ( Wull20lih flows. 

In this paper we determine that for Stokes wave the lowest singularities in of 
both z(w) and n(i(;) are the square-root branch points located periodically at w = 
nX + \Vc — ct, n = 0, ±1, ±2,... (we choose the crests of Stokes wave to be located at 
w = nX — ct) and we determine Vc numerically as a function oi H/X. He re c is the velocity 


of pro pagation of Stokes wave which depends on H. In the previous work iDvachenko et al 
( 20131) we found that as i? —^ Hmax, the branch point approaches real axis with the 
scaling law 


CX {Hmax - Hf 


( 1 . 1 ) 


where 5 = 1.48 ± 0.03. We also provided an accurate estimation of maximum ampli¬ 
tude of the Stokes wave Hmax- Adiabatically slow approach of Stokes wave to its lim¬ 
iting form during wave dynamics is one of the possible routes to wave breaking and 
whitec apping, which are responsib le for significant part of energy dissipation for gravity 
waves ( Zakharov et a/.l[200911200^ . Formation of a close to limiting Stokes wave is also 
considered to be a probable final stage of evolution of a freak (or rogue) waves in the 
ocean resulting in formation of approximate limiting Stokes wave for a limited period of 
time with following wave breaking and disintegration of the wave or whitecapping and 
attenuation of the freak wa ve into wave of regular amplitude ( Rainev fc Longuet-Higgii3 
2006t Zakharov et al.l[2006l) . 

The paper is organized as follows. In Section [5] we introduce the basic equations of 
2D hydrodynamics in conformal conformal variables and reduce these equations to the 
equation for Stokes wave. In Section [3] numerical approaches to simulation of Stokes wave 
are given together with the results of simulations. Also numerical procedures to recover 
the location and type of the branch point are discussed. Section 0] introduces a new 
variable f defining a second conformal transformation which maps one spatial period of 
Stokes wave into the entire real line. Then the Fade approximation of Stokes wave is 
found in complex C plane. T he efficient Alpert-Greengard-Hagstrom (AGH) algorithm 
( Aloert et al. 2000t Lau 2004 1 is used to obtain the Fade approximation. That algorithm 
allows to avoid the appearance of artificial zeros and poles of Fade approximation and 
achieves a spectral accuracy. The convergence of the Fade approximation to the branch 
cut singularity is established which allows to recover the jump at branch cut. Section [S] 
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relates jump at the branch cut in ^ variable to the sum of periodically located branch 
cuts in w complex plane. It is shown how to use the series expansion of the jump along 
branch cuts near branch points to recover the square-root singularity at the branch 
point. It is demonstrated that there are no more singularities in the finite complex plane 
beyond one branch point w = ivc per period. In Section [6] the main results of the paper 
are discussed. Appendix provides a derivation of basic hydrodynamic equations in 
conformal variables. Appendix [B] gives a short description of AGH algorithm adapted for 
Stokes wave. Appendix [C] describes a notation used for the tables of Fade approximants 
for Stokes wave and gives samples of such tables. A full set of tables is provided in the 
electronic attachment. These tables reproduce the Stokes wave with the relative accuracy 
of at least 10“^®. 


2. Basic equations 

In physical coordinates (x, y) a velocity v of 2D potential flow of inviscid incompressible 
fluid is determined by a velocity potential $(x,y,t) as v = V<i). The incompressibility 
condition V • v = 0 results in the Laplace equation 

= 0 ( 2 . 1 ) 


inside fluid —oo < y < r]{x,t). To obtain the closed set of equations we add the de¬ 
caying boundary condition at large depth <l?(x, y, t)|y_>._oo = 0, the kinematic boundary 
condition 


dr] _ f dr] 9 $ 
dt \ dx dx 


and the dynamic boundary condition 



y=r](x,t) 






+ gr] = 0 


y=r](x,t) 


at the free surface 


y = 

We define the boundary value of the velocity potential as 


( 2 . 2 ) 


(2.3) 


(2.4) 


(2.5) 

Consider a time-dependent conformal transformation 


z = z(w,t), w = u + iv 


( 2 . 6 ) 


which maps a half-strip ——c»<r;^0 of complex plane w into a region 
— -I ^ X < —oo < y ^ r]{x, t) of complex physical plane z = x-\-iy each time t such 

that the line —^ u < ^, z; = 0 is mapped into a line of free surface x + ir]{x,t) with 
—-I < X < -I and 


X 





(2.7) 


Also w = —ioo maps into z = —ioo. Here the flow is assumed to be periodic in the 
horizontal direction with the period A both in w and z variables. Conditions suggest 
to separate z{w,t) into a periodic part z{w,t) and a non-periodic part w as follows 


z{w, t) = w -\- z{w, t), or x{w, t) = u + x(w, t), y{w, t) = v + y(w, t), 


( 2 . 8 ) 
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where 


z{w + A) = z{w), X ( ±— ) = 0. 


(2.9) 


Equations (j2.8l) and (1211) extend conformal transformation (I2H) into C . Also x{u,t) 
and y{u,t) form a parametric representation (over the parameter u) of the free surface 
elevation ()2.4I) . 

The idea of using time-dependent confor mal transformation fo r unsteady fluid flow 
was exp l oited by se veral authors including Ovsyannikov (Il973 ). Meison et al. ( 1981 ). 


^nveeijj 1991 j_ 19^ ), ari d Dyachenko et a,l. ( 1996 fJ l: Zakharov et al\ 1 2002b ). We follow 
Dyachenko et a,li (jlQQBoh to recast the system (I2i])-(I2H) into the equivalent form for 


x{u,t), y(u,t) and ip{u,t) at the real line w = u oi the complex plane w using the 
conformal transformation (1^ (see Appendix lAlfor more details). A kinematic boundary 
condition (12.21) is reduced to 

( 2 . 10 ) 
( 2 . 11 ) 


ytXu - xtyu + HiIju = 0 

and the dynamic boundary condition (12.3|) is given by 

f’tyu - f’uyt + gyvu = -H {iptXu - f’uXt + gyxu) ■ 

where 


Hf{u) = ip.v. 


p-\-oo 


fin’) 


du' 


( 2 . 12 ) 


u — u 

is the Hilbert transform with p.v. meaning a Cauchy principal value of integral. Period¬ 
icity of f{u) allows to reduce the integration in the Hilbert transform as follows 

r-A/2 T /.A/2 


Hfiu) = ;^ XI P- 

7T ^' 


fin') 


I -du' = --p.v. , / / s 

l-x/ 2 n' -u + nX A tan (Tt^i^) 


du'. (2.13) 


The equivalence of equations (12.101) and (12.111) to equations (I2.1() - (l2.3p uses the ana- 
lyticity of z(w) and n(ri;) in C“, where 


H = $ -f i© 


(2.14) 


is the complex velocity potential. Here 0 is the stream function defined by 0^ = —4>y 
and 0y = to satisfy Cauchy-Riemann conditions for analyticity of Il{z,t) in z plane. 
The conformal transformation (12.6p ensures that 

0„ = 0„ = (2.15) 

in w plane. The periodicity of the flow implies the condition 

n(u; -I- A, t) = n(u;, t) (2-16) 

together with equation (|2.9I) . We also assumed in equations (j2.10p and (12.111) that 

A/2 A/2 

J ri{x,t)dx = j y{u,t)xuiu,t)du = 0, (2-17) 


-A/2 


-A/2 


meaning that the elevation of free surface of unperturbed fluid is set to zero. The equation 
(|2.17l) is valid at all times and reflects a conservation of the total mass of fluid. 

Both equations (12.101) and (12.111) are defined on the real line w = u. The Hilbert 
operator H transforms into the multiplication operator 

iHf)k = isign(fc)/fe, 


(2.18) 
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for the Fourier coefficients (harmonics) fk, 

xji 


fk 


1 

A 


f{u )exp 



du, 


of the periodic function f{u) 


-A/2 


= f{u + A) represented through the Fourier series 
27T^ 


fiu)= fk exp (iku 


k— — c 


A 


(2.19) 


( 2 . 20 ) 


Here sign(/c) = — 1, 0, 1 for fc < 0, k = 0 and fc > 0, respectively. 

The Fourier series (I2.20p allows to rewrite f{u) = /(w)|^=o as follows 

fiu) = f^{u) + f~{u) + /o, 

where 

OO 

f^{w) = ^/fcexp 
k^l 

is the analytical function in C^, 

-1 

/"M = X! 

k— — oo 




( 2 . 21 ) 

( 2 . 22 ) 


(2.23) 


is the analytical function in C“ and /o = const is the zero harmonic of Fourier series 
(12.201) . In other words, equation (12.211) decompose f{u) into the sum of functions f^{u) 
and f~{u) which are analytically continued from the real line w = u into C’*' and C~, 
respectively. Equations (I2.18|) . (12.211) . (12.221) and (12.231) imply that 

Hf(u)=i[f+(u)-f-(u)]. (2.24) 


If function f(w) is analytic in C“ then f(w) is analytic in C“'" as follows from equations 
(j2.21|) - (j2.23|) . where bar mean complex conjugation, w = u — iv. Then the function f{u), 
u £ M has analytic continuation into because at the real line w = w. Using equations 
dll]) and (12.141) we obtain that = i[n(M,t) +n(u,t)]. It means that after solving 

equations (12.101) and (12.111) one can recover the complex potential H from the analytical 
continuation of 


n(M, t) = 2p'if{u, t) 


(2.25) 


into C . Here 

P=\{l + iH) (2.26) 

is the projector operator, Pf = f~ + ^, into a function which has analytical continuation 
from the real line w = u into C“, as follows from equation (12.241) . Note that without loss 
of generality we assumed the vanishing zero Fourier harmonic, Hq = 0, for n(M,t). 

Also 


H^f = -f (2.27) 

for the function f{u) defined by (12.211) provided the additional restriction that /o = 0 
holds. In other words, the Hilbert transformation is invertible on the class of functions 
represented by their Fourier series provided zeroth Fourier harmonic /o vanishes. If /o 0 
then the identity (12.271) is replaced by 


H^f = -if-fo). 


( 2 . 28 ) 
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The analyticity of z{w) in C 
equations (I2.24|I . (I2.28I1 . that at the real line w = u the following relations hold 


implies, together with x = ^(z + z), y = ■^{z — z) and 


y — yo = Hx and x — xq = —Hy. 


(2.29) 


Here xg and yg are zero Fourier harmonics of x{u,t) and y{u,t), respectively. Note that 
the additio n of zero harmoni c s xg an d into equation (12.291) is the modification compare 
with Refs. lOvachenko et all I 1996q 1: IZakharov et al\ (|2002qI 1 . These Refs, were focused 
on the decaying boundary conditions r]{x,t) —>■ 0 and ip{x,t) —>■ 0 for \x\ —>■ oo which 
imply, together with the condition (I2.17P in the limit A —>■ oo, that io = yo = 0. However, 
generally ig and yg might be nonzero for the periodic solutions with a finite A considered 
in this paper. 

Equations (12.291) imply that it is enough to find either y(u, t) or x{u, t) then the second 
of them is recovered by these explicit expressions. Taking derivative of equations (I2.29|) 
with respect to u results in the similar relations 


Vi 


= Hxu and 


= -Hyu 


(2.30) 


2.1. Progressive waves 

Stokes wave corresponds to a solution of the system (I2.10p and (12.lip in the traveling 
wave form 


1 p{u, t) = ip{u — ct), z(u, t) = z(u — ct), 


(2.31) 


where both ip and z are the periodic functions of u — ct. Here c is the phase velocity of 
Stokes wave. We transform into the moving frame of reference, u — ct ^ u, and assume 
that the crest of the Stokes wave is located at u = 0 as in Fig. [T] and A is the spatial 
period in u variable for both ip and 5 in equation (12.311) . We look for the Stokes wave 
which has one crest per period. Higher order progressive wave s are a lso possible which 
have more than one different peak per period Chen fc Saffman ( 1980ll . However here we 
consider only Stokes wave. We recall that the spatial period A is the same in both u and 
X variables as follows from equation (12.91) . In addition, it implies that the phase velocity 
is the same both in u and x variables so that the Stokes wave has the moving surface 
y = ri{x — ct) and the velocity potential ip = ip{x — ct) in physical spatial variables {x, y) 
with the same value of c as in equations (12.311) . The Stokes solution requires y{u) to be 
the even function while x{u) needs to be the odd function which ensures that y = ri{x — ct) 
is the even function. 

It follows from (12.101) and (12.311) (corresponding to substitution ^ —>• — for y and 
Ip and % —>■ — c|^) that Hipu = cy„ and then excluding ip from (12.111) we obtain that 


- c^yu +gyyu +gH[y{l + Xu)] = 0. (2.32) 

We now apply H to (12.321) . use (12.291) to obtain a closed expression for y, and introduce 
the operator k = —duH = V— which results in the following expression 

ioy= (c^fc-i)y-^^ + yfcy^ =0, (2.33) 

where we made all quantities dimensionless by the following scaling transform u —>■ 
uXf2TC, X —>■ x\I2ti, y —?> yA/27T and c is scaled by cg as follows c —)> ccg, where 
Co = \/glkg is the phase speed of linear gravity wave with the wavenumber kg = 27t/A. 
In these scaled units the period of ip and z is 27T. Our new operator k in Fourier space 
acts as multiplication operator, qualitatively similar to H: {kf)k = kfk- 
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3. Numerical simulation of Stokes wave 

We solve (12.331) numerically to find y{u) by two different methods each of them ben¬ 
eficial for different range of the parameter H/X. For both methods y{u) was expanded 
in cosine Fourier series and the operator k was evaluated numerically using Fast Fourier 
Transform (FFT). A uniform grid with M points was used for th e discretization o f 
—7T ^ u < TC. A first method is inspired by a Petviashvili method ( Pet viashvilil 1 19 7^1 



0.1388. The performance of that method for larger values of H/X was limited by the 
decrease of the speed of numerical convergence. 

3.1. Newton CG and Newton CR methods 

For larger H/X we used a second metho d which is t he Newton Conjugate Gradient 
(Newton-CG) method proposed bv lYana ( 20091 2010ll . The idea behind the Newton- 
CG method is simple and aesthetic: first, linearize (12.331) about the current approxi¬ 
mation j/„, assuming that the exact solution can be written as a sum of current ap¬ 
proximation and a correction y = yn + 5yn'- Loy = 0. Then Loyn + LiSyn — 0, where 
Li = —MSyn — (jt{ynSyn) + y-akSyn + Synkyn^ is the linearization of Lq around the 

current approximation and M = —c^k + 1. Second, solve the resulting linear system 
LiSyn = —Loyn for Syn with one of sta ndard numerical method s, in our case it was ei¬ 
ther Conjugat e Gradient (CG) m ethod ( Flestenes fc Stiefel 1952ll or Gonjugate Residual 
(GR) method ( Luenberger 1970ll to obtain next approximation yn+i = yn + 5y. It should 
be noted that monotonic convergence of GG or GR methods is proven only for positive 
definite (semidefinite for GR) operators, while in our case Li is indefinite. Nevertheless, 
both methods were converging (although generally nonmonotonically) to the solutions, 
and convergence was much faster than using GPM. 

Newton-GG/CR methods can be written in either Fourier space, or in physical space. 
We considered both cases, however Newton-GG/CR methods in Fourier space require 
four fast Fourier transforms per CG/CR step, while in physical space it requires at least 
six. For both cases GG and GR we used M as a preconditioner. 

We found that the region of convergence of the Newton-CG/GR methods to nontrivial 
physical solution (12.3311 is relatively (with respect to GPM) narrow and requires an initial 
guess 2/0 to be quite close to the exact solution y. In practice we first run GPM and then 
choose 2/0 for Newton-CG/GR methods as the last available iterate of GPM. 

Because most of our interest was in getting dependence of characteristics of Stokes 
waves on the wave height and the only parameter in equation (j2.33l) is velocity of prop¬ 
agation c, we were calculating waves changing continuously the parameter c and using 
results of computations with previous values of c as initial condition 2/0 for the New¬ 
ton GG/CR iterations. Due to this approach Newton-GG/CR methods converge to the 
nontrivial solution in all cases provided we additionally used the numerical procedure 
described below in Section [321 


3.2. Stokes wave velocity as a function of steepness 

Results of multiple simulations of Stokes wave are shown in Figure |2 where the wave 
velocity c is shown as a function of the dimensionless wave height H /X. This functio n 
is nonmonotonic which is in agreement with previous simulations (^e.g. ISchwartz ( 1974 1: 
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Figure 2. Oscillations of dimensionless velocity of Stokes wave propagation as function of 
steepness obtained from simulations. First three plots from left to right and from top to bottom 
have increasing zoom both in vertical and horizontal axes to focus on oscillations. In the lower 
right corner a plot is scaled by a magnification function fmag{H /\) — l/{30{Hmax — H )+ 1 
to show all simulation data in a single graph while stressing obtained oscillations. 


WilliamsI ( 1981 . 1985 1) and theoretical analysis ( Longuet-Higgins fc Fox ( 1977 . 1978ll ) 


which predicted an infinite number of oscillations. 


We were able to resolve with quadruple precision two oscillations (two maxima and two 
minima) of the propagation velocity as a function of H/X. These oscillations represent 
a challenge for simulation, because propagation velocity is the only parameter in the 
equation (12.331) . Then it is impossible even to go over the first maximum by changing 
continuously velocity of propagation c. This is because after the maximum is reached, we 
has to start decreasing the parameter c. But decreasing of c causes iterations to converge 
to the less steep solution on the left from the maximum (which we already obtained on 
previous steps), instead of steeper solutions to the right from the maximum. 

In order to resolve this issue we used the following approach. Assume that the singu¬ 
larity of z closest to real axis in w complex plane is the branch point 


z ~ ci(w — iVcY 


(3.1) 


for w —> ivc, where ci is the complex constant, Vc > 0 and /3 are real constants. By the 
periodicity in it, similar branch points are located at ic = ivc + 2nn, n = ±1,±2 ,... 
(recall that that we already switched to the dimensionless coordinates). We expand z(u) 
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Figure 3. Schematic of contour in the C”*" which allows to determine distance Vc from the 

branch cut to the real axis. 


k—0 

into Fourier series z{u) = 5fcexp(ifcu), where 

— OO 


TT 



— 7t 


are Fourier coefficients and the sum is taken over nonpositive integer values of k which 
ensures both 27T-periodicity of z{u) and analyticity of z{w) in C“. We evaluate 
in the limit fc —> — oo by moving the integration contour from the line —n<u<n 
into C’*' until it hits the lowest branch point (13.111 so it goes around branch point and 
continues straight upwards about both sides of the corresponding branch cut as shown 
by the dashed line in right panel of Fig. [3] Here we assume that branch cut is a straight 
line connecting w = iv^ and +ioo. Then the asymptotic of is given by 


|lfc| oc k^-oo. (3.3) 

This approach was used in our previous work ( Dyachenko et al. 2013h to evaluate dis¬ 
tance Vc of the lowest singularity to the real line. Now our key idea is to push artificially 
the singularity w = ivc toward the real line, thus increasing H/X. It follows from ex¬ 
pression (|3.3I1 that to decrease Vc we can multiply Fourier coefficients of the previously 
obtained Stokes wave solution Zk by exp(afc), where the numerical parameter a is cho¬ 
sen such that 0 < a <C Uc. The result of this multiplication Zk exp(afc) is not a Stokes 
wave solution anymore, but it has higher steepness and not very distinct from the Stokes 
wave solution if a is small enough. After that modification we slightly decrease c from 
previous value and allow iterations of Section [5TT] to converge starting from Zk exp(ak) as 
zero iteration. As we expected, iterations then converge to the solution on the right from 
the maximum. This procedure allowed us to resolve both maxima and one nontrivial 
minimum of c as a function of Ff/A as summarized in Fig. 


3.3. Recovering Vc from the Fourier spectrum of Stokes wave 
To obtain the location of the branch point w = iVc with good precision one has to go 
beyond the leading order asymptotic (13.31) . Next order corrections to the integral (13.21) 
for /3 = 1/2 have the following form 

\zk\ ~ -I- C2\k\~^/'^ P C3\k\~'^/‘^ + C4|fcr®/^ -I- .. .^ k -oo, (3.4) 
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Figure 4. Stokes wave with c = 1.082 (blue dash-dotted line), c = 1.091 (green dashed line) 
and c = 1.0922851405 (dark orange solid line). Corresponding values of H/X are given in the 
legend. Inset shows zoom-in into small values of a:/A near a wave crest. 


where we took into account the expansion of z{w) in half-integer powers (in —n = 
0,1, 2,... beyond the leading order term (EH). 

The numerically obtained spectrum \zk \ of Stokes wave was fitted to the expansion (EH) 
in order to recover Vc and coefficients ci, C 2 , C 3 ,.... The highest accuracy in recovering Vc 
was achieved when the middle of spectrum k ^ kmax/‘2 was used for that fit, there k^ax = 
M/2 is the highest Fourier harmonic used in simulations. kmax/‘2 represent a compromise 
between the highest desired values of k to be as close as possible to asymptotic regime 
fc —>■ 00 and the loss of numerical precision for k —> kmax- We estimated the accuracy of 
the fit by varying values of k used for fitting as well as changing the number of terms 
in the expansion (13.41) . Typically we used 4 terms in (13.41) . Section lT3l discusses the 
comparison of the accuracy of the obtained results with the other methods we used to 
find Vc- 


3.4. Highest wave obtained 


We calculated z{u) with high accuracy for different values oi H/X using computations 
in quad precision (32 digits). Such high precision is necessary to reveal the structure of 
singularities in C+. Fig.|4]shows spatial profiles of Stokes waves for several values oi H/X 
in physical variables {x,y). The Stokes wave quickly approaches the profile of limiting 
wave except a small neighborhood of the c rest. 

As it was shown in the previous paper [Pvachenko et al\ ( 20131) . tails of the spectra 
have asymptotic behavior corresponding to /? = 1/2 in (13.3L which means that we have 
square root br anch cut sing ularit y in C~*~ a ll the time. This is consistent with theoretical 
predictions bv iGrant ( 197^ and Tanveeij ( 199ll) . 


The number of Fourier modes M = 2kmax which we used in Fast Fourier Transform 
(in simulations we expand y{u) in cosine Fourier series to speed up simulations and to 
be memory efficient) for each value H/X increases quickly with the increase of H as Vc 
decreases. E.g., for H/X = 0.0994457 it was more than enough to use 256 modes while 
for the largest wave height 


H'f/ff/^/X = 0.141057778854883208164928602256956 (3.5) 

achieved in simulations we used M = 2^^ « 134 x 10® modes. Due to such high num¬ 
ber of modes, the precision of value (EH decreases by round-off errors in approxi- 
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mately times, i.e. in ^ 4 digits. This extreme case has c = 1.0922851405 and 

Vc = 5.9382441 9892803271779 x 1 0~^. These numbers are the moderate extension of our 
previous work Dyachenko et all ( 2013li by pushing down a lowest value of Vc more than 
twice. Further decrease of the numerical values of Vc can be achieved by both subtracting 
the leading order singularity (IMI from the numerical solution and using the nonuniform 
numerical grid in u which concentrates near u = 0. These numerical approaches are 
however beyond t he scope of this papers. 

Before our workiDvachenko et all ( 2013 1. the nu merical estimates of Hm . nx we re found 
bv lWillian3 ( 1985ll as = 0.141063 and Gandzha fc Lukomskv ( 2007 1 

0.14106 348398. The ot her commonly us ed but less precis e estim ate is /X = 

0.1412 ( Schwart^ll974ll . It was shown in Dyachenko et all ( 2013ll that numerical values 
of Vc in the limit (Hmax — H)/X <C 0 was fitted to the scaling law (11.11) with 


HmaxlX = 0.1410633 ±4-10 


-7 


(3.6) 


The mean-square error for 6 in is ~ 0.04 which offers the exact value S = 3/2 as 
a probable candidate for (HID. The estimate (13.61) suggests that the previous estimate 
is more accurate than Als o is within the ac c uracy of the es¬ 
timate (1^ . Howev er, is ob tained in Ref. ICandzha fc Lukomskv! (120071 1 from the 

Michell’s expansion (lMichelllll893ll of the limiting Stokes wave which ignores the expan¬ 
sion in powers of the irrational number /r = 1.4693457 4... Existenc e of that expansion 
beyond the Stokes power law was established bv iGrant ( 1973 1. Lack of resolving 

that expansion suggests that does not have a well controlled accuracy. In contrast, 

our numerical results are based on Fourier series for non-limiting Stokes wave which 
has well-controlled precision. The difference between (13.61) and the new lower boundary 
estimate (13.51) of the largest is ~ 0.004%. 


4. Fade approximation of Stokes wave 

4.1. Additional conformal transformation and spectral convergence of Fade 

approximation 

To analyze the structure of singularities of Stokes wave we perform an additional con¬ 
formal transformation between the complex plane w = u + w and the complex plane for 
the new variable 

C = tan . (4.1) 

Equation (14.11) maps the strip —Tt < Re(w) < n into the complex C, plane. In particular, 
the line segment — tt < w < tt of the real line w = u maps into the real line (—oo, oo) in 
the complex plane C as shown in Fig. [S] Vertical half-lines w = ±7T -|- iu, 0 < w < 
oo are mapped into a branch cut i < C < ioo- In n similar way, vertical half-lines 
w = ±TT -I- iv, —oo < V < 0 are mapped into a branch cut —ioo < C < i. However, 
27T— periodicity of z{w) (1^ allows to ignore these two branch cuts because z{w) is 
continuous across them. Complex infinities w = ±ioo are mapped into f = ±i. An 
unbounded interval [iUc,ioo), Uc > 0 is mapped into a finite interval [ixc 7 i) with 


Xc = tanh■ 


(4.2) 


The mapping (14.11) is dif ferent from the commonly used (see e.g. Schwartj ( 1974ll : Tanveed 
( 1991 ): IWilliamsI ( 198l[ )'l mapping ( = exp (—iic) (maps the strip —tt ^ Re{w) < n into 
the unit circle). The advantage of using the mapping (14.1|) is the compactness of the 
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Figure 5. Schematic of a second conformal map between the periodic domain in w-plane 
(right panel) into C) = tan(w/2)-plane. Another useful property of this map is representation of 
27t-periodic branch cut from iuc to ioo as a finite length cut from Xc = tan(iUc/2) to i. 



W, number of poles 



Figure 6. (a) An exponential decay of error in Fade approximation of Stokes for 
H/\ = 0.125510247666212033511898125908053 as a function of the number of poles N. (b) 
The density p{x) on the branch cut sampled at (^k = iXfc = itanh (^), k — 1,,,N, obtained 
from the Fade approximation of Stokes wave from(a) with N = 29 (green stars). Blue dotted 
line is the estimated profile of p(x) for ths same Stokes wave in the continuous limit of —>■ oo. 


interval (ixc,i) as mapped from the infinite interval (iuc,ic»). In contrast, the mapping 
to the circle leaves the interval (iVc,ioo) infinite in f plane. 

We use Alpert-Greengard-Hagstrom (AGH) algorithm ( Aloert et aZ.ll2000l: Laul2004il to 
approximate the Stokes wave z(C) at the real line Re(C) = C by a set a poles in the complex 
f plane. Approximation by a set of poles is a particular case of Fade approximation by 
rational functions where P{Q and Q(C) are polynomials. Zeros of Q{C) give the 


location of poles. Looking at complex values of f in the rational function provides 
the analytical continuation of z((') into the complex f plane. Usually Fade approximation 
is numerically unstable because of the pairs of spurious zeros and poles appear in finite 
precision arithmetics. These doublets correspond to positions of zeros of P{f) and Q{() 
which are nearly cancel each other. In our practical realizations, AGH algorithm avoids 
the numerical instability of the Fade approximation until the number of poles N increases 
to reach the accuracy corresponding to the round-off error in the numerical approximation 
of z{u). If the analytical continuation of z{u) into w G C has a branch cut, the AGH 
algorithm places poles along the branch cut. AGH algorithm is outlined in Appendix [BJ 
We applied AGH algorithm for z(C) at the real line Re(C) = (, where z{() is obtained 
from simulations described in Section [3l Increasing N we observed the exponential con- 
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Figure 7. The density p(x) for three different Stokes waves in log-log scale. A straight dashed 
line shows scaling law which corresponds to the limiting Stokes wave. Insert shows p(x) in 
linear scale for the same three Stokes wave which are visually almost indistinguishable. 


vergence of Fade approximation z{C)pade to z{C,) as 


err^c, oc e 


-p{vc)N 


(4.3) 


where erroo = max |5(^) — z{Qpade\ is the error in infinity (maximum) norm. An 

— 00<(^<i30 

example of the exponential convergence is shown in Figure [6^ for a particular Stokes 
wave. Here pivc) is the function of Vc but is independent on N. We found that with high 
precision 

p(vc) oc uy®. (4.4) 

AGH algorithm is looking for poles in the entire complex plane (. All the encountered 
poles for Stokes wave were found on the interval of imaginary axis along the interval 
[iXc,i), where Xc is determined numerically as in Section |3l 
Equations (14.31) and (14.41) demonstrate excellent performance of Fade approximation. 
E.g., decreasing Vc by six order required in our simulations only 10-fold increase of N 
as detailed in Appendix C. It suggests that numerical method which solves Stokes wave 
equation (j2.33|) directly in terms of Fade approximants might be superior to Fourier 
methods including numerical approaches mentioned in Section [3.31 This topic is however 
beyond the scope of this paper. 

It is rather straightforward to distinguish in AGH algorithm poles from branch cuts. If 
both poles and branch cuts would be present in z(C) then increasing N one observes that 
some poles of Fade approximation are not moving and their complex residues remain 
approximately the same. These correspond to poles of 5(^). Such behavior occurs for 
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test problems when we artificially added extra poles to z{(). Other poles of Fade ap¬ 
proximation are moving with the increase of N and their complex residues are changing. 
These poles mark the spatial location of branch cuts of 5(C). The density of poles along 
each branch cut is increasing with the increase of TV. If the jump of 5(C) at branch cut 
is continuous along it then we expect to see the convergence of density of poles with the 
increase of N. All this is valid until err^ decreases down to the level of round-off error 
at which 5(C) was determined. Further increase of N would result in the appearance of 
spurious poles at random positions of C plane with the magnitudes of complex residues 
at the level of round off error (^ 10“^^ for z{C,) found with quad precision in Section[3]). 

Using 5(C) obtained by the method of Section [31 we found a single branch cut [ixc,i) 
but no poles in Stokes wave. It means that in complex w plane we have one branch cut 
per spatial period 27t located at {2nn -|- ivc, 2nn + ioo), n G N. 

We parametrize that branch cut as follows 





pjx'W 

C-ix' 


(4.5) 


where p(x) is the density along branch cut. That density is related to the jump of 5(C) 
at branch cut as explained in Section 15.11 The constant yo is determined by the value 
of 5(C)|c=oo = z(w)\w=Ti- This constant has a zero imaginary part, Im(?/o) = 0, because 
i('u;)|u,=n /2 = 0 as given by the equation (12.911 . 

The Fade approximation represents equation (14.511 as follows 





Pix')dx' 

C-ix' 


lyo 


N 


^iC-ix. 


(4.6) 


where the numerical values of the pole positions Xn and the complex residues 7 „ (n = 
1,..., N) are obtained from AGH algorithm. 


4.2. Recovering jump along branch cut 

We recover p{x) from equation (14.611 as follows. Assume that we approximate the integral 
in equation (14.511 by the trapezoidal rule 


1 



Xc 


Pix')dx' 
C - ix' 


X2 - Xi Pi 

2 C - iXi 


N-l 




Xn-\-l Xn—1 Pn 

2 C - iXn 


Xn — Xn-1 Pn 

2 C - iXJV ’ 


(4.7) 


A comparison of equations (14.611 and (14.711 suggests the approximation pn,N of the density 
piXn) on the discrete grid Xn, n = 1,..., iV as follows 

p{Xn) ze pn,N = -- for n = 2,..., iV - 1, (4.8a) 

Xn+1 Xn—1 

p{Xi) - Pi,N = p{xn) ze pn,n = -- ■ (4.8&) 

X2 - Xi XN - XN-l 

A convergence of pn,N to the continuous limit p(xn) as N increases is quadratic with 
the error scaling oc for x away from boundaries X = Xc and x = 1- Near these 
boundaries we cannot apply the trapezoidal rule and have to resort to less accurate 
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estimates given by the equation (j4.86p . Figure [8] demonstrates this oc convergence of 
the Fade approximation to the continuous limit. Figure |6 }d shows the particular example 
of p{x) (shown by solid line solid line) compared with pn,N (shown by stars) for N = 29. 
We believe that the convergence of to the continuous value p{x) SuS N ^ oo and 
absence of other poles outside of [iycj i] provide a numerical proof that the only singularity 
of z(C) are the branch points ixc and i connected by the branch cut ( G [ixc,i]- 

At X = Xc the function p(x) has a square root singularity as given below by equation 
(j4.13ll . This singularity additionally reduces the accuracy of the approximation (|4.8bll for 
pi^N which is based on Taylor series. To significantly improve numerical accuracy of 
we assume that p has the following square root dependence in the vicinity of Xc'- 

Papproxix) ~ A\/X Xc- (4'9) 


Here the values of the parameters A and Xc are determined from two interior points 
{X 2 i P 2 ,n) and (xa, P 3 ,n) found via the trapezoid rule (14.8al) . We assume that Papprox(X 2 ) = 
P 2 ,N and Papproxixa) = P 3 ,N which gives that 


f P3,N - P2,N \ P3 ,nX2 - pInXS 

V X3-X2 J P3,N-P2,N 


(4.10) 


Using equations (14.91) and (14.101) for x = Xi we obtain the numerically accurate approx¬ 
imation that 


Pl,N 


( (X3 - Xi)pI,n - (X2 - Xi)pI,n 
\ X3 - X2 


1/2 


(4.11) 


where p 2 ,Ar and p^^^ are given by equation (I4.8al) . 

At X = 1 the function p{x) also has singularity and respectively numerical value of 
Pn,n from (14.861) is not very accurate. To improve that accuracy we use that p{l) = 1 
as found in Part II. Then using the trapezoidal rule we obtain much more accurate 
expression that 


p{Xn) — Pn,n = — 


2xn 


Xn-1 


(4.12) 


Figure 0 shows the density p(x) for three different Stokes waves in log-log scaling. It 
is also seen that inside the branch cut and for small Xc ^ 1, the density p{x) scales as 
X^/^ which corresponds to the limiting Stokes wave. A deviation from that scaling occurs 

near x = Xc and x = 1- __ 

Classical Markov’s theorem ( Markofllll895 ) proves pointwise convergence of the diag¬ 
onal Fade approximants [N/N]f of the function / of the type (14.51) with p(x) ^ 0 in 
the limit fV —>■ oo for () G C \ [ixc,i]. Here the diagonal Fade approximation [N/N]f 
of the function / means that both polynomials P(C) and Q{C) has the same order N 
which is natural for the discretization (14.61) . More general Fade approximants of the func¬ 
tion / are [N/M]f, wher e N and M are the orde r s of t he polynomials P(C) and Q{C), 
respectively. Theorem of de Montessus de Ballord ( 1902ll ensures pointwise convergence 
of [N/M]f —>■ / for —>■ oo with fixed M in the disk |C| < i? if / the meromorphic 
function in that disk with exactly M poles (counted according to their multiplicity). 
However, the diagonal Fade approximations of the meromorphic function / generally 
fails t o provide u n iform convergence wi th the known counterexamples given by Buslaev! 
( 2001 ): Lubinsk-^l ( 20031) . Nuttalll ( 19701 ) showed that instead the diagonal Fade approx¬ 
imants of meromorphic function for N ^ oo have a weaker convergence in logarithmic 
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N 

Figure 8. Error between Fade approximation with N poles and the continuous limit for p{x) for 
Stokes wave with H/X — 0.125510247666212033511898125908053. It is seen that the error oc 
for large N. To calculate that error we use a spline interpolation for with N = Nmax = 33 to 
construct the approximation of the continuous limit of the square of the density, pcontinuousix)- 

After that the error is defined as err = - plontinuousiXn)f') /{N -2) for each 

N, where p„,jv is given by equation (ITSI) . 


capacity which allows the lack of pointwise convergence along exceptional sets. Gonchar 
I Gonchailll973L 1975) extended Markov’s theorem on the pointwise convergence of the 
diagonal Fade approximants to the functions f + r, where / is the function of the type 
(j4.5ll with p(x) > 0 almost everywhere in % € [xc, 1] and r is the meromorphic function 
away from branch cut and has no poles at branch cut. Convergence in logarithmic capac¬ 
ity of the diagonal Fade approximants of the analytic function /(C) with a finite number 
of branch points (this is a more general type than the type (14.51) because thes e branc h 
points can be located away from a single line) was proven bv IStahl ( IDSbol fft. 1997 ). 
That convergence o ccurs away from c ertain sets of C (in some cases these sets are simple 
arcs). See also Ref. Aptekarev et al\ ( 201ll) for the recent review. All these results were 
obtained for Fade approximants based on the Taylor series at a single point in C. Thus 
these results do not directly apply to AGH algorithm which is based on least squares 
approximation at mu ltiple points of C. AGH algorithm is also distin ct from multipoint 
Fade approximation ( G. A. Baker fc Graves-Morrisl[l9^ Safi 1972 ). where the Taylor 
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series is interpolated at multiple points in contrast to least squares in A GH algorithm. 


Fade approximants were also constructed based on least squares in Ref. iGonnet et al 


(j201lh were it was conjectured that least squares-type algorithms might ensure pointwise 
convergence. That conjecture is consistent with our simulations. 


4.3. Finding a numerical value of a location of branch point C = ixc 

There are different ways to find the location of branch point ^ = iXc from simulations. 
First way is based on the decay of Fourier spectrum of z{C,) for n » 1 and is described 
in Section [3731 Second way is to find p{xn), n = 1,... ,N and then determine the point 
p{Xc) = 0 by the polynomial extrapolation of p{xn)- First and second ways provide 
comparable numerical accuracy in our simulations (typically the relative error in Xc is 

- 10-4). 

We found however, that better accuracy is achieved in the third way as follows. Con¬ 
sider the formal series 

OO 

Zser = X! (4.13) 

j=o 

in the neighborhood of the branch point C, = iXc- The term ieb"/4 in front of the coef¬ 
ficients Oj is chosen for convenience to ensure that coefficients Oj take real values. The 
radius of convergence of that series is 2xc as discussed in Part II. Taking M = 10 — 20 
terms in that series one can use the nonlinear fit to determine the unknowns Xc and Oj. 
Typically we use Nj = 30 — 40 points («„, z(u„)) such that all values are inside the 
disk of convergence |u„ — ixc\ < ‘^Xc of the series (14.131) . Here values of 5(u„) are taken 
from simulations of Section [3] with being the numerical grid points closest to m = 0. 
The accuracy of the nonlinear fit is typically ^ lO-^^ as estimated by varying M and 
Nm- In Part II we provide much more accurate way of calculating Xc which is based on 
the compatibility of the series (14.131) with the equation (I2.33|) of Stokes wave. In con¬ 
trast, the above three methods use numerical values of («„, z(un)) obtained as described 
in Section [3] and do not use the equation ()2.33(1 directly. 


5. Stokes wave as an integral over jump at branch cut and the 
expansion of density p near a branch point 


5.1. Jump at branch cut 


Sokhotski-Plemelj theorem (see e. 
plied to the equation (14.51) gives 


g. iGakhovI ( 1966ll : IPolvanin fc Manzhir^ ( 2008 11 ap- 


z{ix ± 0) = iyo + p.v. f ± np{x), Xc < X < 1- (5.1) 

J i(x-x) 

Xc 

Thus the jump of z(() at branch cut is —2np(x) for crossing branch cut at C = iy in 
counterclockwise direction. 


5.2. Stokes wave as the sum of contribution from branch cuts in w complex plane 

Gonsider a representation of Stokes wave by the density p along branch cuts in complex 
plane w. Because of the 27T-periodicity in u direction we write z{w) as the integral over 
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periodically located branch cuts, 


OO 

5H=^i+/ E ( 

n= —OO ^ 


w + 27xn — iv' b + 2nn — iv' 


p{v')dv', 


(5.2) 


where zi is the complex constant, Vc is related to Xc by ( 021 ) , a summation over n ensure 
the periodicity of z{w) along u and we replaced p{x) by p{v') to distinguish it from p(x) 

in (14.51) . Also we introduced the additional term-which is intended to 

b + 2nn - iv' 

ensure the convergence of the integral. The constant b can be chosen at our convenience. 
A change of that constant results in the change of Zi. 

The sum in (15.2p is then calculated using of the identity 


E 


1 

n + a 


7T cot na 


giving 


z{w) = Zi + ■ 


cot 


W — IV 


— cot 


b-iv' 


p{v')dv'. 


(5.3) 


Taking the limit Im{w) —> —c» we obtain from equations (15.31) and (14.51) that 


z{u — i(X)) ~ + 2 / ' * ~ 


p{v')dv' = iyo + 


P(x)dx 

-i-ix' 


(5.4) 


We set 


and 


X =tanh — 
2 


P(x) = p(2arctanhx). 


(5.5) 


(5.6) 


We also require that Zi = iyo then we find from the equations (15.41) . (15.51) and (15.6p 
that 


b — 71. 


Using the trigonometric identity 


cot {a — b) = 


1 + tan a tan b 


tan a — tan b 

one obtains from (j5.3p . (15.4L ()5.5p . (j5.6p and (|5.7p that 


CXD 

z{w) = iyo + ^ J 


1 + i tan Y tanh ^ 

-—itanh-^ 1 p(v')dv' 

tan ^ — i tanh ^ 2 


(5.7) 


OO 1 

i?/0 + ^ / ( C-fy “ = iyo + J 


P(x)dx 

C- iX ’ 


(5.8) 


i.e. we recovered the equation (14.51) from the equation (15.21) . 
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5.3. Expansion of p{x) in powers of f — ixc 
Assume that we have the branch cut (ixcA) ^iC) in th® complex plane of f and that 
the branch point at C = iXc is of square root type. Then we expand p(x) in the following 
series 


p(x)=E^ 2 n+l(x-Xc)'/^+”. 


(5.9) 


n =0 


Note that that adding terms of integer powers of (C — iXc) into the equation (15.91) is not 
allowed because it would produce logarithmic singularity at ^ = ixc throug h the equatio n 
(14.51) wh i ch is incompatible with the Stokes wave as was shown in Refs. iGrant (1973); 
Tanveeil ( 199lh . 


Integrating over x in ()4.5I1 using ()5.9p gives 

/(C) = h ( 2i^l - Xc - 2iV Xc + iC arctan 


\/l - Xc 


VXc + iC. 


+ h ( g\/l - Xc(i - 4ixc + 3C) + 2i(xc + iC)^^^ arctan 


\/l - Xc 


Vxc + K. 

+ bs - Xc (3i - lliXc + 23iXc + 5C - 35xcC - 15iC^) 


— 2 i(xc + iC)^^^ arctan 


VI - Xc 


VXc + iC. 


br(. ■.) 


(5.10) 


A series expansion of (IS.lOp at C = iXc and comparison with the series (14.131) result in 
the relations 


by+i = (-ir+^ay+i, j =0,1,2,... (5.11) 

Note that the expansion (15.101) provides the relations for bn with only odd values of n. 
This is because the series dSH) is convergent only inside its disk of convergence, x~Xc < 'r, 
where r is the radius of convergence. It will be shown in Part II that r = 2xc for Xc < 1 /3. 
The explicit expression for p(x) is unknown for Xc + ^ < X < 1 while /o(x) still contributes 
to the terms a2j{C — iXc)V J = 0,1, 2,... in the series (14.131) . 

Thus the expansion (15.9p together with the relations (15.111) provides a convenient tool 
to work with p(x) near to x = Xc- 


5.4. Absence of singularities in branch cut beyond the branch points f = iXc cmd f = i. 

A priori one can not exclude existence of singularities inside the branch cut f G [ixc, i] 
beyond branch points f = lYr and d = i at its ends. Existence of such singularities 
were conjectured in Refs. [Grant ( 1973 ): ISchwart j ( 1974ll . To address that possibility we 
subtracted the expansion (15.101) from the numerical solution of z{C,) for Stokes wave. We 
obtained both z{f) and recovered p(x) through AGH algorithm using variable precision 
arithmetics with ^ 200 digits to achieve a high precision in that subtraction. Typically 
we used Stokes wave of the moderate nonlinearity with Xc ^ 10“^ to operate with 
the moderate number of required Fourier harmonics. After that the numerical values of 
^ 1 : b^, & 5 ,... were recovered from fitting of p{x) to the expansion (15.91) near x = Xc- 
Typically we truncated the expansion (15. 9p to the first 3 terms bi, 63 , 65 which results 
in the truncated function f (Cf)truncated in the expansion (I5.10L Also Xc was obtained by 
the procedures described in Section IT^ Alternative way to recover bi, 63 , 65 is through 
using the expansion (I3a was also used but generally gives lower precision. 

Next step was to take mth derivative of ziC,) — f{C)truncated over f numerically and 
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obtain the Fade approximation for the resulting expression [z(^) — f (Otruncated]^^'^ re¬ 
sulting in new density p{x)- If singularity would be present inside the branch cut 
then it would correspond to singularity in p{x)- However, we did not find any sign of 
such singularities at least for moderate order of derivative m = 1,2,3. It suggests that 
f = iXc and f = i are the only si ngularitie s in co mplex f plain. This conclusion is in 
agreement with the results of both iTanveen (1199111 and Part II obtained by alternative 
methods. 


6. Conclusion 

In this paper we found numerically the Stokes solutions of the primordial Euler equa¬ 
tions with free surface for large range of wave heights, including the approach to the 
limiting Stokes wave. The limiting Stokes wave emerges as the singularity reaches the 
fluid surface. We found from our high precision simulations (between 32 and more than 
200 digits) the Fade approximation of branch cut singularity of Stokes wave. We provided 
the tables of the Fade approximants for a wide range of Stokes wave steepness. These 
tables allow to recover Stokes wave with the minimum accuracy 10“^®. We show that 
these Fade approximants quickly converge to the jump at branch cut as the number of 
poles N increases with the scaling law (ig.(l4Hl). We use the series expansion of the 
jump along branch cuts in half integer powers to recover the square-root singularity at 
the branch point. We found that there are no more singularities in the finite complex 
plane beyond one branch point per period. Following Part II is devoted to the analysis 
of the structure and location of branch points in infinite set of sheets of Riemann surface 
beyond the physical sheet of Riemann surface considered here. 

The authors would like to thank Prof. S. Lau for the introduction to AGH method of 
Fade approximation and sharing his computer codes which were used at the initial stage 
of research. Also the aut hors thank developers of FFTW (|Frigo fc Johnson! 120051) and 
the whole GNU project (IGNU Project! 11984-20121) for developing, and supporting this 
useful and free software. The work of S.D. and A.K. was partially supported by the U.S. 
National Science Foundation (grant no. OGE 1131791). The work of A.K. and P.L. on the 
Fade approximation was supported by Russian Science Foundation grant 14-22-00259. 


Appendix A. Derivation of dynamical equations 

In this Appendix we adapt the work of iDvachenko et~^ (1996a) to the case of the 
periodic boundary conditions deriving the basic dynamical equations (12.101) and (12.lip 
for 2D ideal hydrodynamics with free surface in conformal variables. We use similar 


notations to 


Dyachenko et al\ (11996011 . 


Dyachenko et al\ 1 1996 all and provide steps of the derivation skipped in 


A.l. Hamiltonian after conformal map 

It was shown by IZakharov ( 1968 1. that the potential flow of an ideal fluid with free 
surface is the canonical Hamiltonian system with canonical variables rj (1^ and Ip (USD. 
Ganonical Hamiltonian equations 


dr] 5H dfj SH 

dt Sip ’ dt Sr] 


(Al) 


are equivalent to the boundary conditions (12.21) and (12.31) . Here H is the Hamiltonian 
which coincides with the total energy (the sum of the kinetic energy T and the potential 
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Figure 9. A schematic of one period of a wave with a counterclockwise contour of integration 

for application of Green’s theorem. 


energy U) per spatial period of wave A, 


H = T+U = 


A/2 n A/2 

^ / dx / (V$)^dy+| J ri'^dx, 

-A/2 -oo -\/2 


and without loss of generality the fluid density is set to one. One has to express the 
kinetic energy 

A/2 r) 

T=^ J dx y'(Vd>)2dy (A2) 

-A/2 -oo 

through the canonical variables 77 and tp which generally requires to solve the Laplace 
equation (EU with the boundary conditions (12.21) . (12.3L (12.51) and $(a;, y, t)|i /->—00 = 0 
in the region ^ x < ^, —00 < y ^ r]{x,t). That region is schematically shown in 
Figure [9l Using relations 


V . (4.V<I.) = (V*)^ + , (V4.)“ = A (*|i) + A (*|) , 


which are valid for the harmonic function $ (EH) and applying Green’s theorem to the 
equation (lA 21) one obtains that 


2T = 


A/2 T) 



-A/ 2-00 C 


^9$ *9<i) 

dx + $^dj/ 
oy ox 


(A3) 


Here C is a positively (counterclockwise) oriented contour along the boundary of the 
periodic domain occupied by fluid shown in Figure |9l A sum of integrals along left and 
right hand sides (vertical segments) of the contour vanishes due to periodicity. Integral 
along lower part of the contour (horizontal segment) is zero due to the boundary condition 
on potential $j,=_oo = 0. Notice that in the case of finite depth fluid with a rigid flat 
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bottom y = —h, the integral along similar lower segment y = —h is also zero because the 
boundary condition at the hnite depth bottom is ^y\y=-h = 0 (zero vertical velocity at 
the bottom) and dy = 0. Then the equation (lA 3|) is reduced to the following line integral 


r:=-\j2 


2T = 


x—Xj2, y—r){x,t) 


—dx 
dy 




(A 4) 


We use the time-dependent conformal transformation (EH),(1121) to relate partial deriva¬ 
tives in cc, ?/ and u,v as follows 


which implies that 



d<^ 


9$ 

du 

OX 

+ 

"q 

dy 





dv 

ox 

+ 

dy 


^uXu 


^vVu 

dx 

rr2 

+ 

9 ’ 

yi 


^uVu 

+ 


dy 

x^ 

+ 

2 ’ 
yi 


(A 5) 
(A 6) 


where we also used Cauchy-Riemann equations = y^ and Xy = —for the conformal 
map z{w). 

Substituting (|A 51) and (IA6I) into (IA4|) . using relations dx = Xydu and dy = y„dit on 
the line w = u one obtains that 


2T = 


-A/2 A/2 

J = J 




(A 7) 


A/2 -A/2 

Here we took into account the orientation of the contour and conditions EH on confor¬ 
mal transformation. 

Sokhotski-Plemelj theorem (|5.ip (seee.g. Gakhov ( 19661) : IPolvanin fc Manzhirovl l 20081) 1 
allows to express a real part of the function which is analytic in the lower (upper) half 
plane through the imaginary part (and vice versa) at the real line u = w using the 
Hilbert transformation (|2.12[) . For a conformal transformation z{w, t) = x{w, t) -\-iy{w, t) 
such relations are given by (12.291) . A complex velocity potential n(z,t) = $ -|- i0 is the 
analytic function in the fluid domain —oo < y ^ y(a:, t), where 0 is the stream function. 
The conformal transformation z = z(ic, t) (12.61) ensures that H remains analytic function 
after transforming from z to w variable with the lower half plane C“ being the domain 
of analyticity in w. Similar to equation (I2.29|) . real and imaginary parts of H are related 
at the real line u = w through the Hilbert transformation as follows 


0 = $ = -HQ, 


(AS) 


where we assumed the decaying boundary condition = 0|t,=-oo = 0. Here we 

abuse notation and use the same H and $ both for independent variables w and z : 
^{w,t) = ^{z,t) and fl{w,t) = n(z,t), i.e. we omit tilde. 

Also the analyticity of H implies that the velocity potential <f> is the harmonic function 
satisfying the Laplace equation EH both in X, y variables and similarly 

V^<i>('u, v,t) = 0 
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in variables u and v. Using Cauchy-Riemann equations (12.1511 and the relations (jA 81) one 
obtains that 

= (A 9) 

Substituting (IA9|) into (lA 71) we express the kinetic energy in terms of canonical variable 
ijj as follows 

A/2 A/2 

2T = [ 4)„$|^^pdu = - f (A 10) 

-A/2 -A/2 

Here we used the definition (1^ which in w plane turns into ip(u,t) = = 0,t) 

as follows from the mapping of the fluid surface into the real line u = 0. Then the 
Hamiltonian in terms of variables on the surface takes the following form 

A/2 A/2 

i? = -i f tljHtljudu+^ ( y'^Xudu. (All) 

-A/2 -A/2 


A.2. Least action principle in conformal variables 

We use the constrained Lagrangian formulation to obtain the dynamical equations in 
conformal variables at fluid surface. A time dependence of the map (I2.6I) implies that 
we have to ensure the analyticity of that map through the appropriate constraint. We 
discuss the Lagrangian dynamics first and add the corresponding constraint later in this 
Section. Equations dAH) realize extremum of an action 


S = 



(A 12) 


with the Lagrangian 



U—da; - H. 
at 


(A 13) 


The first term here has to be converted from the integral over x into u variable. Consider 
mapping (x, t) —>■ {u, r), which is the change of parametrization of the surface under the 
conformal map. Here t = t. Transformation u = u{x, t) is the inverse to the conformal 
map X = x(u,t). The fluid surface r]{x,t) after transformation corresponds to y(u,T). 
We express drj/dt by the chain rule as follows 


dr] 


dy dr dy du 
dr dt du dt 


dy dy du 
dr du dt' 


(A 14) 


To find du/dt here we express full differentials of x and t through u and r as follows 


f dx \ Xu Xr \ f du \ _ -j- / du \ 

[ dt )-[ tu u )[ dr )=’^[ dr )■ 

Taking into account that t = t, one obtains the Jacobian matrix 

j _ f \ 

0 1 )■ 


(A 15) 


(A 16) 
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Inverse procedure for full differentials of u and r through x and t yields that 

f \ — f du/dx du/dt \ / da; \ _ i / da; \ 

V dr j “ V dr/dt ) \ dt ) ~ \ ' 

Comparing entries of the matrix in (lA 171) with inverse of (lA 161) . one gets 

du Xt 

dt Xu 
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(A 17) 


(A 18) 


Substituting (lA 181) into (lA 141) yields 

dr] Xt 


(A 19) 


We use the Lagrangian (lA 131) to substitute it into the action (lA 121) . Consider the first 
term in the action, 

S = J jLdt = j Jfj^dxdt + ... (A20) 

and perform a change of variables in the integral as dxdt = det(J)dudT = a;ududT. 
Together with the expression (jA 19|) it results in 


^gdxdi = 



Xududr 


fpiyrXu - yuXT)dudT. (A 21) 


Using (lA 111) . (lA 131) . (IA21I) and adding the analyticity constraint (12.291) ensuring that 
y — yo = Hx and taking into account that t = t as well we obtain a new constrained 
Lagrangian 


A/2 


A/2 


A/2 


L= / V.HVvd.-| / !,=x.d. 


(A 22) 


-\/2 


-A/2 


-A/2 


A/2 


+ / {y-yo- Hx)fdu, 

-A/2 

where / is the Lagrange multiplier for the analyticity constraint. 


A. 3. Variations of action 

We now obtain the dynamical equations from the Hamilton’s least action principle. Van¬ 
ishing of variational derivative SS/6ijj = 0 of the action (jA 12p with the Lagrangian (lA 22p 
over potential ip on the surface yields the following expression 

ytXu - yuXt + HiIju=0- (A 23) 

This equation is nothing else but kinematic boundary condition (1^ after the conformal 
map into w plane. 

Two conditions 6S/dx = 0 and 6S/Sy = 0 result in equations 

yu'ipt - ytipu + gyyu = Hf, (A 24) 

-Xuipt + xtipu - gyxu = /, (A 25) 

which are turned into a single equation by excluding the Lagrange multiplier / giving 
yuf’t - ytf’u + H{xuipt - xtipu) + g[yyu + H{yxu)\ = 0. (A 26) 


Equations (lA 231) and (IA26I) recover the implicit dynamical equations (12.1011 and (12.1111 . 
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A.4. Zeroth harmonic in implicit dynamical equations (I2.10|) and (j2.1ip and 

conservation of momentum 

Consider Fourier transformations of the surface elevation y{u, t) and the velocity potential 
on surface if with respect to conformal coordinate u, 


y{u,t) = yoi^t) + 

k^O 

if{u,f) = ipo{t) + ipkit)e^^'^. 
fc/O 

Here zeroth harmonics yo(t) and ifoit) are written separately and are given by 

A/2 A/2 


(A 27) 


yoit) = 


y(u,t)du, tfoit) = 


ip{u, t)du. 


-A/2 -A/2 

One can rewrite equation (lA 261) in the following form 

d 

yui’t - yttpu = -H{xuift - xti>u + yxu) - -^-^y^- 

2 au 


(A 28) 


(A 29) 


A zeroth Fourier harmonic of the right hand side (r.h.s.) of equation (lA 291) vanishes 
because the term in parenthesis is multiplied by H which removes any zeroth harmonic 
and the remaining term is the partial derivative over u. Respectively, the zeroth harmonic 
of the left hand side (l.h.s.) of equation (lA 291) must vanish. Integrating that l.h.s. to 
obtain the zeroth harmonic, using equation (lA 291) and integrating by parts over u one 
obtains that 

A/2 A/2 A/2 

/ {vuf’t - yti^u)du = - / {yulft + yut'f’)du = 


1 

A 


A 


A dt 


iiyudu = d, (A 30) 


-A/2 


-\/2 


-A/2 


A/2 

where we used a periodicity of ip and y in u. Thus f ipyudu is the integral of motion. 

-A/2 

To find a physical meaning of that integral we note that natural candidates for conserved 
quantities are the components of the total momentum of fluid along x and y directions. 
Taking into account that fluid density is one, we obtain x component of momentum 
as an integral of the horizontal velocity inside fluid, which gives 


A/2 77 ( 21 ,i) —A/2 

Px= J j = J = J 

-A/2 -00 C A/2 


, dy 


A/2 


dx = — 


y=r](x,t) 


ipyudu. (A 31) 


-A/2 


Here we applied Green’s theorem to positively oriented contour C shown in Figure [9l 
Due to periodicity of functions and decaying boundary condition $(a;, y, t)|y_i._oo = 0, 
only integral along the surface is nonzero. Comparison of equations (lA 301) and (IA3ip 
shows that consistency of equation (lA 291) is ensured by the conservation of the horizontal 
component P^ of the total momentum of the fluid. 

Applying the Hilbert transformation H to equation (IA29I) and using the identity (12.281) 
one obtains that 




xtipu + yxu - go = H 



yttpn 


g d . 
2^^ 


(A 32) 


where go is the zeroth Fourier harmonic of Xuipt — xtipu + yXu- To find go we proceed 
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similar to equations (jA 30|) and (|A31(1 to find that 

A/2 A/2 


A/2 A/2 A/2 

J Mt-xti^u + yxu)du=~ j i^Xudu+j J yxudu, (A33) 

-A/2 -A/2 -A/2 

\/2 

where J ipXudu is the integral of motion corresponding to the conservation of the 

-A/2 

vertical component Py of the total momentum of fluid, 


Py = 


A/2 A/2 

/ dx / *.d„ = /-*dx= j d.x,d„. 


(A 34) 


-A/2 -00 C -A/2 

Then equations (lA 33|) and (lA 341) imply that qq is the integral of motion given by 

A/2 A/2 


qo = j J yxudu=j J ri{x,t)dx 


(A 35) 


-A/2 


-A/2 


and representing a conservation of the total mass of fluid. Also according to equation 
(|2.17|) . we set go = 0 in this paper. 


Appendix B. Alpert-Greengard-Hagstrom (AGH) Algorithm and 
Stokes Wave 

In this Appendix we describe an efficient algo rithm for Fade app roximation o f the 
functi on on a discrete grid, following original work lAlnert et al\ (2000) and work bv iLad 
( 2004h where more detailed explanation and further development of the algorithm was 
presented. 

Consider 27T-periodic complex-valued function f{u) = z(u) — u — iyo defined on a 
grid with nodes uj € [—tt, 7t]. Values of the function at the grid points are denoted as 
fj = f{uj). We look for an approximation of f{u) in the form of a ratio of two polynomials 
P{u) and Q{u), i.e. the Fade approximation. We briefly describe AGH algorithm in a 
general way with additional comments for our particular case. As it was mentioned in 
Section m we use the second conformal map C = tan(u/2). The introduction of auxiliary 
variable allows to consider the real line C G R as opposed to considering a finite interval 
u G [—7T,Tt], while the infinity along the imaginary axis is mapped into imaginary unit i 
and 27T-periodicity in u direction is ensured. Without loss of generality we assume that 
/(±7t) = 0. In this paper, we take f{u) = z{u) — u — iyo = z(u) — iyo (see equations (|4.5p 
and (14.61) for comparison). Condition /(±7t) = 0 allows to consider P and Q such that 
the degree of polynomials are degQ = 1 -bdegP = N, where the integer N is allowed to 
vary. We are looking for the convergence of the rational approximation to /, 


QiO 


/(O, 


in a sense of solving a minimization problem 


+7T 


mm 
P,Q , 


P{u) 


Q{u) 


- fiu) 


du. 


(Bl) 
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That minimization problem is challenging because Q in the denominator makes (EH 
nonlinear problem. In the transformed variable C the problem (|B 1|) remains nonlinear 
and is reduced to 

+ CX3 


mm 

P,Q 


PiO 


QiO 


-/(C) 


dC 

C + i' 


(B2) 


In AGH algorithm, the complexity of nonlinearity is bypassed by solving instead of 
a sequence of linear least-square problems 


+ 00 


mm 

p(i + l)_Q(i+l) 






/(«) 


du, i = 1, 2,... 


(B3) 


We define an inner product 


+ CX3 

{f,9)i= J f{u)g{u)wt{u)du, 


(B4) 


— OO 
1 


with a weight function Wi{u) = (for ^-plane the formula for the weight is 

modified to be Wi{() = I/(|( 5 d)(^)| 2(^2 _|_ Then the previous least squares problem 
can be rewritten as follows 


where 

Uum = {9,9h 

is the norm. 

As it was shown in Aloert et a.l. ( 2000l l. the solution of the least squares problem 
is equivalent to the solution of 

(-P*+i+g*+V(u),/r„(u)), =0, 

for n = I,..., 2N, with hn{u) defined as follows 

forn = 2,4,6, ...,2Af, 

1 for n = 1, 3, 5,..., 2Af - I, 


(B5) 


(B6) 


which are nothing else but 

f{u), 1, uf{u),u, v^fiu), 






This claim can be proven by variation of mth coefficient of P(u) (for even n = 2m + 2) 
and Q{u) (for odd n = 2m -1-1). We put coefficient at the leading power of Q{u) to be 
equal to one. Thus we need to find 2N coefficients for two polynomials. 

We orthogonalize 2N + 1 functions hn{u) using Gramm-Schmidt orthogonalization 
procedure, 


9n{u) = 


' f{u), forn = l, 

1 - C 2 i/(u), for n = 2, 

min{4,n —1} 

ugn- 2 {u)- X) Cnj9n-jiu), ioT n = 3,... ,2N + 1, 


i=i 


(B7) 
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where real constants Cnj 

are given by 



^nj \ ((ug„_ 2 ,ff„- 3 >i 

for n = 2 and j = 1, 

., for n = 3,..., 2N + 1 and j = 1,.. 

., min{4, n — 1}. 

(B8) 

L {9n — j :9n — j )i 


Then we obtain that 

g2N+i = -P^^+^'> + f{u)Q^^+^\ 


(B9) 


so and are computed from recurrence coefficients Cnj by splitting into even 

and odd-numbered parts. 

For our purposes of finding the jump at branch cut it is convenient to represent a ratio 
of P{u) and Q{u) as a sum of simple poles, 


Pju) ^ ^ 7n 
Qiu) 


(BIO) 


In order to do that we compute zeros Xn, n = 1,2,...,TV of Q(u) using Newton’s 
iterations. At each step one zero Xn is found by Newton’s iterations. After that we 
remove that zero from Q{u) by division on [u — Xn) and proceed to the next step for the 
modified Q etc. After that procedure coefficients 7 „ are given by the following expression 


P{Xn) 


Q'iXn)' 


(Bll) 


Derivative Q'(u) are obtained from previous recurrence relation for gn{u) by differentia¬ 
tion. 


Appendix C. Tables of Stokes Waves 

Using the Fade approximation, introduced in Section lU one can approximate Stokes 
wave for each value of the scaled height P[/X as a sum of poles 


N 

z{w) ~ Zpade[u) = W + lyo + 

n—1 


In 

td,n{w/2) - ix„ 


(Cl) 


Here N is the number of poles in the Fade approximation. Using AGH algorithm (see 
Appendix IbJ we found that all poles for all values of H/\ are located on the imaginary 
axis. 

We provide Tables [T]l4] for four particular cases of Stokes waves with wave heights 
ranging from HjL ~ 0.031791 to PdfL ~ 0.141058. Complete library of computed 
wav es can be accessed throu gh the electronic attachments as well as through the web 
link Dyachenko et al\ ( 2015ll . These data of Fade approximation allow to recover the 
Stokes wave with the relative accuracy of at least 10“^® (for the vast majority of cases 
the actual accuracy is higher by several orders of magnitude). First and second columns of 
both Tables and electronic files represent values of Xn and 7 „, respectively. Additionally 
a third column in electronic files provides the values oi pn,N, n = 1, 2,..., A calculated 
from data of the first two columns using equations (I4.8al) . (14.111) and (14.121) . 

We used three quantities to characterize the accuracy of our numerical Stokes wave 
solution and its Fade approximation. First quantity is the residue 


M 


1/2 


R{y) = N ( XI l^oy(uj)P 

D=i 
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of equation (j2.33|) . R{y) characterizes convergence of our iteration algorithm described 
in Section o to the Stokes wave. Here M = 2kmax is the number of grid points uj used 
in the discretization of z{u). Second quantity is the relative error of Fade approximation 


^'^'^pade 


/ M 

S ~ Zpade{Uj)\'^ 




M 


E 

\ 


1/2 


of our numerical solution z(uj). Third quantity is the amplitude of the highest Fourier 
harmonics lifcmaxl used in FFT. 

We balanced these three quantities in our simulation to achieve the most efficient and 
reliable approximants of Stokes waves. Typically we chose k^ax large enough such that 
\zk I < to ensure that our discretization error is below 10“^®. Here the 

factor characterize the accumulation of round-off error in FFTs. A convergence 

of numerical iterations down to R{y) — 10“^® was found to be sufficient to achieve the 
desired accuracy of solution in 10“^®. After that we used AGH algorithm with N large 
enough to make sure that ervpade is below 10“^® by several orders of magnitude. 

The second and third rows in electronic . dat-files provide the additional information 
extracted from simulations which include the number of points of the numerical grid 
M = 2kmax, the residual R{y), the Stokes wave height y^ sX x = ±7T, the amplitude 
of the highest Fourier harmonics |.Zfcmax|j the Fade error err pads, the scaled Stokes wave 
height H/X and the Stokes wave velocity c. Values of H/X are also encoded in the names 
of .dat-files. Also the file summary.txt provides a summary of the results from all .dat- 
files. 
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i.96041092606335083862992746108661e-01 
i.78972925087544517288851755005498e-01 
i.49569603918982534434611327659588e-01 
i.l0406118678767801011884022998371e-01 
;.64694768844023775849318292632706e-01 
;.15784392644967788264370902239031e-01 
■.66774518804464747111211133286936e-01 
'.20283901206785220595281417979068e-01 
1.78365420413127130751057141705342e-01 
1.42527484790841967661585477231496e-01 
1.13814667765069562012985702014292e-01 
..929083157742315710201020784180286-01 
..802208826392953721044026130453576-01 


7.86955267798815779896940975384730e-03 

1.589382086499895499702201560075586-02 

2.09270477914666067462762444957813e-02 

2.307728553092327749219270329785936-02 

2.27821823395535616299282467454014e-02 

2.06781569334250898322582073617042e-02 

1.74505018291868390667201976549862e-02 

1.371915186205552703870528555878686-02 

9.979034893934588248110470184706106-03 

6.58739944600110438403394267024205e-03 

3.78022620936003042036390618792491e-03 

1.70068438916655758680561664969251e-03 

4.279395361919985118980051772408956-04 


Table 1. Data for Fade approximation of tho wav6 with volocity 
c = 1.005, th6 st66pn6ss H/\ = 0.031791185830078550217424174610939, and 

yo = —0.094819818875344225940453182945545. Paramotors of simulations and Fade ap¬ 
proximation are M = 16384, R{y) — 3.64 x 10“®®, ervpade — 4.65 x 10“®®, and tho smallost 


Fourior harmonic had valu6 \zk„ 


~ 1.00 X 10" 
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k Xk 

1 9.95104877443162988543285604210300e-01 

2 9.74036453796113160099814623502153e-01 

3 9.37570348097817646693771937993553e-01 

4 8.88487837568099082583936075213862e-01 

5 8.30226745019764721513358604279726e-01 

6 7.66333992991599455804305516283483e-01 

7 7.00044654389407424698796692053902e-01 

8 6.34030062620564487061621925382604e-01 

9 5.70306175038391450182479612388228e-01 

10 5.10259067527902434999804572718914e-01 

11 4.54736803912594713785686635920502e-01 

12 4.04166063099136862080187337641328e-01 

13 3.58666569676145947416712220270910e-01 

14 3.18149576080072883402812844826148e-01 

15 2.82395821449537219671027831006306e-01 

16 2.51113624911206124013006119628666e-01 

17 2.23980136883660703565243361779229e-01 

18 2.00669410169964201214033391681290e-01 

19 1.80870715726851513117763940777923e-01 

20 1.64299946133899344356900387886524e-01 

21 1.50706307186061197360270221491920e-01 

22 1.39875923302152112463866847539994e-01 

23 1.31633517886184493678824771993353e-01 

24 1.25842975665215746611096122396649e-01 

25 1.22407333541749966033626434463374e-01 


7fc 

9.86344259137750131929660816853428e-03 
2.04507668432927246329100238980833e-02 
2.79392155141519735983551710052556e-02 
3.23466950514503536891484865775305e-02 
3.4014778159425464029321685563311le-02 
3.35031653839781443242921103012960e-02 
3.14489535896823273634285653396772e-02 
2.84500818501058994363745193170182e-02 
2.49972283728050063777939045774063e-02 
2.14511658748252720344184761694938e-02 
1.80506679368664453110659093724408e-02 
1.49353242488736140402002871952586e-02 
1.21719831322217078183473222652247e-02 
9.77859794087821305028150640251530e-03 
7.74304960914297497673324174361282e-03 
6.03673166180751660604462350895498e-03 
4.62368871094252860720889700284368e-03 
3.46637996984113299044046402602488e-03 
2.52906751261377992038440816275726e-03 
1.77962987225519971005973434397249e-03 
1.19038860012342014795358493343618e-03 
7.38354660230910637368329193600582e-04 
4.05164363997062298811235878087573e-04 
1.76877312168739070195683856144061e-04 
4.37430353029634816661994533800963e-05 


Table 2. Data for Fade approximation of the wave with velocity 
c = 1.051, the steepness H/X = 0.10042675172528485854673515635249, and 

yo = —0.25732914098527682158156915646871. Parameters of simulations and Fade ap¬ 
proximation are M = 16384, R{y) — 5.19 x 10“®^, ervpade — 1.69 x 10“®^, and the smallest 
Fourier harmonic had value l^fc^axl — 1-00 ^ 10“®^. 
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k Xk 

1 9.95433825932608550132384034857213e-01 

2 9.75720732729872361661639075445556e-01 

3 9.41458580699036796998398974950063e-01 

4 8.95072079198109934865810125596955e-01 

5 8.39603368860036006278810820983949e-01 

6 7.78246707846125864530269738095348e-01 

7 7.13977463595480003911919153327480e-01 

8 6.49314375127269358025790723963870e-01 

9 5.86215465928989576951120366245581e-01 

10 5.26077937447751872984275728101522e-01 

11 4.69802407403950485965504029996914e-01 

12 4.17886269376110040045868819592273e-01 

13 3.70521444182340763884420027666827e-01 

14 3.27682439755960213510467363606042e-01 

15 2.89198722549172239607721584343698e-01 

16 2.54810466354666195816665370686997e-01 

17 2.24209367334721589574533037466343e-01 

18 1.97067225067065397541246930106441e-01 

19 1.73055091362268451120364081334665e-01 

20 1.51855462026532436104149361379210e-01 

21 1.33169515540203788156642166906591e-01 

22 1.16720933274701291000887725539493e-01 

23 1.02257431405919038178275031736290e-01 

24 8.95508127283858594116407019601820e-02 

25 7.83961028319930614767173088020279e-02 

26 6.86101566854466616210976539289158e-02 

27 6.00299941184905270702652324982135e-02 

28 5.25110331386998700818247974509258e-02 

29 4.59253281231247667007365599867943e-02 

30 4.01598777968588267195919995947710e-02 

31 3.51150396988382900972228874525109e-02 

32 3.07030693071633025856068326424687e-02 

33 2.68467901323811952321185618074578e-02 

34 2.34783937316708681288640025635146e-02 

35 2.05383642178678839820064784655908e-02 

36 1.79745193852614117777103084010265e-02 

37 1.57411593704806090331919997928406e-02 

38 1.37983133839445750806859736203507e-02 

39 1.21110752013028358188007191935072e-02 

40 1.06490185910779922966530559824050e-02 

41 9.38568452652012956402075369290278e-03 

42 8.29813278439420925351354881096343e-03 

43 7.36655130339489438745177566237321e-03 

44 6.57391741813790189283510916305043e-03 

45 5.90570578080460414802523913159400e-03 

46 5.34963842759485046283496838777743e-03 

47 4.89547304507788970804144262729396e-03 

48 4.53482604688247134485690783763108e-03 

49 4.26102758786433892841273275905455e-03 

50 4.06900612679011263031777402788905e-03 

51 3.95520060861743412410650474994064e-03 


7fc 

9.06513968659594263994154983859022e-03 
1.88104172479867625607988377665013e-02 
2.58024268240062770737774505111002e-02 
3.00760315302764871790369706833779e-02 
3.19167851791100543244698807489031e-02 
3.17900776761535529310595879964877e-02 
3.02335467307095385878230877915997e-02 
2.77625384158746106085149228773401e-02 
2.48107842884599707048520815911264e-02 
2.17068351144352237760178371853073e-02 
1.86762941525117699745376670289500e-02 
1.58579969144004349124961258254260e-02 
1.33248510736740059585812992162263e-02 
1.11037283880977232781708686935009e-02 
9.19185285105063318623096617861130e-03 
7.56907825149316224628365974410675e-03 
6.20644453419544059242384682030945e-03 
5.07177162315234989673990118919696e-03 
4.13307962854491967554778057762908e-03 
3.36050786167674783997799826083550e-03 
2.72724894316322072325277221844458e-03 
2.20986855953690886419588367706623e-03 
1.78826390264321999398279771128206e-03 
1.44542657328550726857080252749622e-03 
1.16711442195804019558665126539427e-03 
9.41495432894175937842874927421205e-04 
7.58799737152752549854592147665995e-04 
6.10998699362372783997691390696487e-04 
4.91519483298208817983041233383980e-04 
3.94997269772740759120747730149340e-04 
3.17063815302569437242780768838413e-04 
2.54169263655189066241987434556100e-04 
2.03433380026332773186591672337583e-04 
1.62522237707382077096302354822532e-04 
1.29546576264629788075276452826328e-04 
1.02978399018687148172885868196623e-04 
8.15827862763857084235776103594481e-05 
6.43623147351034758281773572995289e-05 
5.05118646282732166085507134151195e-05 
3.93819507325944359184298112513636e-05 
3.04490256349179913250724576956079e-05 
2.32914730544885187107604807304053e-05 
1.75702379530898436181417357806391e-05 
1.30132325361793585312305562950935e-05 
9.40281755981180973510534726928533e-06 
6.56579130822932954380944114441017e-06 
4.36542870899759471322771653831421e-06 
2.69520431445801848314129686051970e-06 
1.47390871215592819437170943388163e-06 
6.41931271126488105657466176599801e-07 
1.58535598642490907942230532511450e-07 


Table 3. Data for Fade approximation of the wave with velocity 
c = 1.0929, the steepness H/X = 0.13825830866311310404416736817381, and 

yo = —0.2915339172431288292999965032009. Parameters of simulations and Fade ap¬ 
proximation are M = 65536, R{y) — 2.59 x 10“®^, ervpade — 1.01 x 10“^^, and the smallest 
Fourier harmonic had value \zkmax \ — 5-00 10“®®. 
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k Xk 

1 9.93398643583003025153415435504531e-01 

2 9.65060482058669453480870980252106e-01 

3 9.16714789789161849202896809849705e-01 

4 8.53133231672953078297597841861823e-01 

5 7.79861194757267505357891461462460e-01 

6 7.02176870736814955937036169062899e-01 

7 6.24431910013545641111933110534714e-01 

8 5.49795332847499273525022520666518e-01 

9 4.80296232546974823630387324524545e-01 

10 4.17024645857800791824977075950470e-01 

11 3.60378305224974615221811203424356e-01 

12 3.10289842646985418345605736936048e-01 

13 2.66407648286693861932841258994721e-01 

14 2.28226886286517332438415939592608e-01 

15 1.95177843545134733789119725180096e-01 

16 1.66681805023593130823224346797399e-01 

17 1.42183993584453030868317439497003e-01 

18 1.21171172484979583446779037786714e-01 

19 1.03179452511723992079582272565214e-01 

20 8.77961228595568606401648108610398e-02 

21 7.46580315447597465030574978595406e-02 

22 6.3448130644394071950746312766621le-02 

23 5.38911872290903942028496746096495e-02 

24 4.57492587887819300483008251206986e-02 

25 3.88172752742240105088258746631741e-02 

26 3.29189096590597195144947958198357e-02 

27 2.79028212558501660025568098114343e-02 

28 2.36392981742008042750164949870111e-02 

29 2.00172924085211268884110103807234e-02 

30 1.69418235041002160106748854919734e-02 

31 1.43317184931477833930656499987995e-02 

32 1.21176530092033068365142201489852e-02 

33 1.02404588173178159766150762312400e-02 

34 8.64966499042977177325090332415729e-03 

35 7.30224274879133050003140619052446e-03 

36 6.16152705009717845830463720403643e-03 

37 5.19629108350495067218636168760723e-03 

38 4.37995272297094390520221610919454e-03 

39 3.68989465744363535075551189083628e-03 

40 3.10688231013106004394611827804953e-03 

41 2.61456578557885722141496603598660e-03 

42 2.19905395393225001547396749458147e-03 

43 1.84855041619759608964615915855160e-03 

44 1.55304251639257745230445095694940e-03 

45 1.30403580024761683344961061119599e-03 

46 1.09432738764763509113618247317394e-03 

47 9.17812647645417749658202151298786e-04 

48 7.69320359699922614599799922262156e-04 

49 6.44472229750674596538655661872796e-04 

50 5.39563219268408452082296554789566e-04 

51 4.51459652187210048535645581625337e-04 

52 3.77512500082907183102947952504772e-04 

53 3.15483620098192985576130479232556e-04 

54 2.63483041310020855659728424130368e-04 

55 2.19915670752623096428691515132722e-04 


7fc 

1.28741415762741679346145664996829e-02 

2.60267721567181859054226183119919e-02 

3.43697679627825352269776411757291e-02 

3.81940153445550033591099941584920e-02 

3.83579555843352602595202964927154e-02 

3.59774931820106663192752740784546e-02 

3.21282856711997437340001627375168e-02 

2.76691730121067711930378384900232e-02 

2.31906560517515347893552651709825e-02 

1.90429324287558759377242844714520e-02 

1.53959152585563261217878597243142e-02 

1.23006351687210032636876990145911e-02 

9.73852232172353639036354649203610e-03 

7.65583101348162009153247587551195e-03 

5.98535170044127311114464478095979e-03 

4.65887488461708640776107914945203e-03 

3.61357577121517253016005783742974e-03 

2.79470260489814969664959294926652e-03 

2.15617836676118577227243251241331e-03 

1.66012679310648786275121618740345e-03 

1.27592107645668605269203574613739e-03 

9.79089087450639485345366472519338e-04 

7.50248095624531518737269429111365e-04 

5.74148509163471182463355721181721e-04 

4.38854383779301056123215099079188e-04 

3.35061475940689863018834966768472e-04 

2.55540871176311150876065830758264e-04 

1.94691329540109355414952726522057e-04 

1.48182722053010736525844218770677e-04 

1.12674155188771601572543993641392e-04 

8.55924207962470484084537963575900e-05 

6.49586463790799447674548874436684e-05 

4.92531498561398293539370127249816e-05 

3.73103914016284265439688122490895e-05 

2.82375225733516039477667066021446e-05 

2.13513661747602139469739101816769e-05 

1.61297466879220448079982238240246e-05 

1.21739653477311336991048098474804e-05 

9.17991094830884893227377698101798e-06 

6.91584935818872179036239759592572e-06 

5.20536938490202499152488268561648e-06 

3.91430351823738560080383189178292e-06 

2.94070847979427725012698794615673e-06 

2.20719808861136536998364371692094e-06 

1.65508232204597074734186170013816e-06 

1.23989163441506187230543611758337e-06 

9.27962520159008244607471765496508e-07 

6.93835226000040035242065190155287e-07 

5.18272340903010891482796683060429e-07 

3.86751562364387865498882802928547e-07 

2.88320243535854856966320396641389e-07 

2.14725698524118689197774299181126e-07 

1.59755495420848083506690922212921e-07 

1.18737501808356300532888071897965e-07 

8.81613507227033650558982811933814e-08 
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k Xk 

56 1.83436026617218827592565025305802e-04 

57 1.52909808789403664334781521734162e-04 

58 1.27381290492618707318705788922866e-04 

59 1.06045663430404223545931459157107e-04 

60 8.82255960623097328938047290772380e-05 

61 7.33513735582244917463595672161340e-05 

62 6.09440811020089107359438382849930e-05 

63 5.06013718010847795103344003355303e-05 

64 4.19854284842886757235756347132190e-05 

65 3.48127867887149603598695647961378e-05 

66 2.88457365344826502539670719441546e-05 

67 2.38850607015635551072991988129380e-05 

68 1.97639074377385155071521574022480e-05 

69 1.63426213283602794258356970211462e-05 

70 1.35043863909173156500338815494453e-05 

71 1.11515555862156824939722316951789e-05 

72 9.20256064170640198246888785687659e-06 

73 7.58931213612359462965716406185546e-06 

74 6.25501349165233025416447266351498e-06 

75 5.15232409179580442223607363260130e-06 

76 4.24181649775016522927709678588896e-06 

77 3.49068106812952109834665722037669e-06 

78 2.87163870904281978123381120499882e-06 

79 2.36202927278644614134938592885535e-06 

80 1.94304885285494843909466767533574e-06 

81 1.59911321065017380641753434790850e-06 

82 1.31732715240169671433477796333342e-06 

83 1.08704211469651017023210815498690e-06 

84 8.99487103028405911079240475052776e-07 

85 7.47460944734835494603150605681252e-07 

86 6.25076012803183072949714996778214e-07 

87 5.27545188407840404034972402364704e-07 

88 4.51005279547543069512840045259275e-07 

89 3.92371966727686847652422460137122e-07 

90 3.49224206935226307895172607412510e-07 

91 3.19719913320901893897552159916508e-07 

92 3.02547325679113057732806265306883e-07 


7fc 

6.53921058513165146118265854353209e-08 
4.84538716133603951388247789431925e-08 
3.58664158964311161459118236867907e-08 
2.65219345830887370991540659742850e-08 
1.95921833026360368361816879242692e-08 
1.44585603393432845280965460758372e-08 
1.06595193757093831390156057934750e-08 
7.85105510061063345075831379965022e-09 
5.77704094551584320668305953789402e-09 
4.24698455656734771613508619208412e-09 
3.11936272695846692633031613888051e-09 
2.28914896008025877293425583671707e-09 
1.67848842336096800302384900691726e-09 
1.22973329115351088076908027814457e-09 
9.00246494394826687766241103823572e-10 
6.58529018137743573464498106258884e-10 
4.81336779425064219659348751354589e-10 
3.51536622960117080648972753857693e-10 
2.56513789011984140244514391180857e-10 
1.86990411531843046778670738725091e-10 
1.36150027070583877598207373400005e-10 
9.89896125536513948820894793870582e-ll 
7.18405497951763524448150948318652e-ll 
5.20148665829541462756399079612847e-ll 
3.75443719949043320491081089015406e-ll 
2.69886771259819847047462342227678e-ll 
1.92941913934894250545144244154976e-ll 
1.36906826573455708305877340415078e-ll 
9.61543204537063425390253166033986e-12 
6.65765514720154023428202200838344e-12 
4.51781968589147435681661138798625e-12 
2.97790805796735353056158063859767e-12 
1.87970510537824614741588711602232e-12 
1.10896242549626643385666064787054e-12 
5.83901681870167816802893811429244e-13 
2.47055052842115889755229649961340e-13 
5.99062680825025191284334866467265e-14 


Table 4. Data for Fade approximation of the wave with velocity 
c = 1.0922851405, the steepness H/\ = 0.14105777885488320816492860225696, and 

yo = —0.289784811618456872977429611644. Parameters of simulations and Fade approxi¬ 
mation are M = 134217728, R{y) — 6.14 x 10“^^, ervpade — 5.43 x 10“^^, and the smallest 
Fourier harmonic had value l^fc^axl — 3-00 ^ 10“®^. 
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